#include "graphsupporttest/TestPrimalDualFactor3Solver.h"

#include <sstream>
#include <numeric>

#include <boost/random.hpp>

#include "boosttestsupport.h"

#include "accumulsupport/AccumulSupport.h"

#include "graphsupport/solve/facilitylocation/PrimalDualSolution.h"
#include "graphsupport/solve/facilitylocation/IFacilityLocationSolver.h"
#include "graphsupport/solve/facilitylocation/PrimalDualFacilitiesLocationSolver.h"
#include "graphsupport/solve/facilitylocation/BinaryRunFacilityLocationSolver.h"
#include "graphsupport/graph/operation/MakeSymmOperation.h"
#include "graphsupport/graph/operation/TriangleInequalityCloseOperation.h"
#include "graphsupport/graph/operation/GenBipartLinkNear.h"
#include "graphsupport/graph/operation/CompoundGraphOperation.h"
#include "graphsupport/graph/operation/AddChpExpWeightsOperation.h"
#include "graphsupport/graph/operation/GeneratePricesOperation.h"
#include "graphsupport/graph/generate/DefaultGraphGenerator.h"

void TestPrimDualSolvContext()
{
	FacilitiesGraph::FacilityPriceType EPS = 1e-5;

	int cityCount = 30;
	int fcltCount = 10;
	PrimDualSolvContext context( cityCount, fcltCount );
	std::fill( context.fgraph.facilitiesPrices.begin(), context.fgraph.facilitiesPrices.end(), 0 );

//	BOOST_CHECK_EQUAL( cityCount, context.cityPays4total.size() );
	BOOST_CHECK_EQUAL( cityCount, context.cityPays4fclt.size() );
	BOOST_CHECK_EQUAL( fcltCount, context.cityPays4fclt[0].size() );

	Graph::AdjMatrixType adjMatr = context.fgraph.getAdjacentMatrixArr();

	for ( int i = 0; i < cityCount; i++ )
	{
		for ( int j = 0; j < fcltCount; j++ )
		{
			Graph::AdjCoefType roadPrice = context.getRoadPrice( i, j );
			BOOST_CHECK_LE( context.getCityPays4road( i, j ), roadPrice - EPS );
		}
	}

	for ( int j = 0; j < fcltCount; j++ )
	{
		FacilitiesGraph::FacilityPriceType totalPayed = 0;
		for ( int i = 0; i < cityCount; i++ )
		{
			totalPayed += context.cityPays4fclt[i][j];
		}
		BOOST_CHECK_LE( context.fgraph.facilitiesPrices[j], totalPayed + EPS );
	}
}

void TestFindWitnessForVertex()
{
	try {
	FacilitiesGraph::FacilityPriceType EPS = 1e-5;

	int citiesCount = 70;
	int facilitiesCount = 30;
	PrimDualSolvContext scontext( citiesCount, facilitiesCount );

	Graph::AdjMatrixType adjMatr = scontext.fgraph.getAdjacentMatrixArr();

	adjMatr[0][citiesCount + 0] = 123;
	scontext.fgraph.facilitiesPrices[0] = 567;

	adjMatr[0][citiesCount + 1] = 234;
	scontext.fgraph.facilitiesPrices[1] = 678;

	adjMatr[0][citiesCount + 12] = 345;
	scontext.fgraph.facilitiesPrices[12] = 123 + 567 - 345 - 1;

	adjMatr[0][citiesCount + 23] = 456;
	scontext.fgraph.facilitiesPrices[23] = 789;

	int witness = PrimalDualSolution::createWitnessForCity( scontext, 0 );
	CRITICAL_CHECK( witness >= 0 );
	CRITICAL_CHECK( witness < facilitiesCount );
	BOOST_CHECK_EQUAL( 12, witness );
	for ( int i = 0; i < facilitiesCount; i++ )
	{
		BOOST_CHECK_LE( scontext.getPayedForFacility( i ), scontext.fgraph.facilitiesPrices[i] + EPS );
	}
	BOOST_CHECK_CLOSE( scontext.getCityPays4road( 0, witness ), scontext.fgraph.getFSEdge( 0, witness ), EPS * 100 );
	BOOST_CHECK_GT( scontext.cityPays4fclt[0][witness], 0 );

	adjMatr[1][citiesCount + 5] = 123;
	scontext.fgraph.facilitiesPrices[5] = 567;

	adjMatr[1][citiesCount + 7] = 234;
	scontext.fgraph.facilitiesPrices[7] = 678;

	adjMatr[1][citiesCount + 12] = 123 + 567 - 1;
//	scontext.fgraph.facilityPrices[12] = was set before

	adjMatr[1][citiesCount + 25] = 456;

	scontext.fgraph.facilitiesPrices[25] = 789;
	witness = PrimalDualSolution::createWitnessForCity( scontext, 1 );
	CRITICAL_CHECK( witness >= 0 );
	CRITICAL_CHECK( witness < facilitiesCount );
	BOOST_CHECK_EQUAL( 12, witness );
	for ( int i = 0; i < facilitiesCount; i++ )
	{
		BOOST_CHECK_LE( scontext.getPayedForFacility( i ), scontext.fgraph.facilitiesPrices[i] + EPS );
	}
	BOOST_CHECK_CLOSE( scontext.getCityPays4road( 1, witness ), scontext.fgraph.getFSEdge( 1, witness ), EPS * 100 );
	BOOST_CHECK_CLOSE( scontext.cityPays4fclt[1][witness], 0, EPS * 100 );
	} catch ( SolutionProcessException & exc )
	{
		std::cout << "\n\n" << exc.message << "\n\n";
		BOOST_CHECK( false );
	}
}

void Test_FindAllWitnesses_RandomCases()
{
	FacilitiesGraph::FacilityPriceType EPS = 1e-5;

	int maxPrice = 300;
	int fracPow = ( int )( 1e2 ); // 2 digits

	int testCaseCount = 5;
	int citiesCount = 300;
	int facilitiesCount = 50;
	double noEdgeProbPerc = 30; // 30%

	boost::variate_generator< boost::mt19937, boost::uniform_int<> > sampler( boost::mt19937( ( unsigned ) 12345 ), boost::uniform_int<>( 0, 1 << 30 ) );

	for ( int tc = 0; tc < testCaseCount; tc++ )
	{
		PrimDualSolvContext scontext( citiesCount, facilitiesCount );
		Graph::AdjMatrixType adjMatr = scontext.fgraph.getAdjacentMatrixArr();

		for ( int i = 0; i < citiesCount; i++ )
		{
			for ( int j = 0; j < facilitiesCount; j++ )
			{
				double noedge = sampler() % 100000; // 100%
				if ( noedge < noEdgeProbPerc * 1000 ) continue;
				FacilitiesGraph::FacilityPriceType price = sampler() % ( maxPrice * fracPow );
				price /= fracPow;
				adjMatr[i][j + citiesCount] = price;
			}
		}
		for ( int j = 0; j < facilitiesCount; j++ )
		{
			FacilitiesGraph::FacilityPriceType price = sampler() % ( maxPrice * fracPow );
			price /= fracPow;
			scontext.fgraph.facilitiesPrices[j] = price;
		}

		bool created = PrimalDualSolution::createWitnessForAllCities( scontext );
		CRITICAL_CHECK( created );

		for ( int i = 0; i < citiesCount; i++ )
		{
			int witness = PrimalDualSolution::findWitnessForCity( scontext, i );
			BOOST_CHECK_GT( witness, -1 );
			if ( witness == -1 ) break;
			for ( int j = 0; j < facilitiesCount; j++ )
			{
				FacilitiesGraph::FacilityPriceType cityPays4road = scontext.getCityPays4road( i, j );
				FacilitiesGraph::FacilityPriceType edgePrice = scontext.fgraph.getFSEdge( i, j );
				BOOST_CHECK_LE( cityPays4road, edgePrice + EPS );
				BOOST_CHECK_CLOSE( scontext.fgraph.facilitiesPrices[witness], scontext.getPayedForFacility( witness ), EPS * 100 );
			}
		}
		for ( int j = 0; j < facilitiesCount; j++ )
		{
			FacilitiesGraph::FacilityPriceType payedForFacility = scontext.getPayedForFacility( j );
			FacilitiesGraph::FacilityPriceType facilityPrice = scontext.fgraph.facilitiesPrices[j];
			BOOST_CHECK_LE( payedForFacility, facilityPrice + EPS );
		}
	}
}

void TestComputeSpecialEdgesSquareSubgraph()
{
	int citiesCount = 100;
	int facilitiesCount = 100;
	PrimDualSolvContext scontext( citiesCount, facilitiesCount );

	scontext.fgraph.setFSEdge( 0, 10, 123 );
	scontext.fgraph.setFSEdge( 1, 10, 234 );
	scontext.fgraph.setFSEdge( 1, 20, 345 );
	scontext.fgraph.setFSEdge( 2, 20, 456 );
	scontext.fgraph.setFSEdge( 2, 30, 567 );
	scontext.fgraph.setFSEdge( 3, 30, 678 );
	scontext.fgraph.setFSEdge( 3, 40, 789 );

	scontext.fgraph.setFSEdge( 0, 40, 998 );

	MakeSymmOperation::makeSymm( scontext.fgraph );

	scontext.fgraph.facilitiesPrices[10] += 998 - 123;
	scontext.setCityPays4total( 0, 998 );

	scontext.fgraph.facilitiesPrices[10] += 346 - 234;
	scontext.fgraph.facilitiesPrices[20] += 346 - 345;
	scontext.setCityPays4total( 1, 346 );

	scontext.fgraph.facilitiesPrices[20] += 568 - 456;
	scontext.fgraph.facilitiesPrices[30] += 568 - 567;
	scontext.setCityPays4total( 2, 568 );

	scontext.fgraph.facilitiesPrices[30] += 789+1 - 678;
	scontext.fgraph.facilitiesPrices[40] += 789+1 - 789;
	scontext.setCityPays4total( 3, 789+1 );

	PrimalDualSolution::computeSpecialEdgesSquareSubgraph( scontext );

	BOOST_CHECK( scontext.specialSubgraph.existsEdge( 10, 20 ) );
	BOOST_CHECK( scontext.specialSubgraph.existsEdge( 20, 10 ) );
	BOOST_CHECK( scontext.specialSubgraph.existsEdge( 20, 30 ) );
	BOOST_CHECK( scontext.specialSubgraph.existsEdge( 30, 20 ) );
	BOOST_CHECK( scontext.specialSubgraph.existsEdge( 30, 40 ) );
	BOOST_CHECK( scontext.specialSubgraph.existsEdge( 40, 30 ) );
	BOOST_CHECK( !scontext.specialSubgraph.existsEdge( 10, 30 ) );
	BOOST_CHECK( !scontext.specialSubgraph.existsEdge( 30, 10 ) );
	BOOST_CHECK( !scontext.specialSubgraph.existsEdge( 10, 40 ) );
	BOOST_CHECK( !scontext.specialSubgraph.existsEdge( 40, 10 ) );
}

bool checkIndependent( PrimDualSolvContext & scontext, std::ostream & buff,
	FacilitiesGraph::FacilityPriceType eps = 1e-5 )
{
	bool result = true;
	int facilitiesCount = scontext.fgraph.getFacilitiesCount();

	for ( int fclt = 0; fclt < facilitiesCount; fclt++ )
	{
		if ( scontext.indepFacilities[fclt] )
		{
			if ( scontext.getPayedForFacility( fclt ) < scontext.fgraph.facilitiesPrices[fclt] - eps )
			{
				buff << "\nfacility " << fclt
					<< " is in independent set, but it was not completely payed: price: "
					<< scontext.fgraph.facilitiesPrices[fclt] << " payment: "
					<< scontext.getPayedForFacility( fclt );
				result = false;
			}
		}
	}

	if ( !result )
	{
		return result;
	}

	for ( int i = 0; i < facilitiesCount; i++ )
	{
		if ( scontext.getPayedForFacility( i ) < scontext.fgraph.facilitiesPrices[i] - eps )
		{
			continue;//do not consider such facilities
		}
		for ( int j = 0; j < facilitiesCount; j++ )
		{
			if ( scontext.getPayedForFacility( j ) < scontext.fgraph.facilitiesPrices[j] - eps )
			{
				continue;//do not consider such facilities
			}
			if ( scontext.specialSubgraph.existsEdge( i, j ) && scontext.indepFacilities[i] &&
				scontext.indepFacilities[j] )
			{
				buff << "\nvertexes " << i << ", " << j << " both in set, but exists edge " << i << " -> " << j;
				result = false;
			}
		}
	}

	for ( int i = 0; i < facilitiesCount; i++ )
	{
		if ( scontext.indepFacilities[i] )
		{
			continue;
		}
		if ( scontext.getPayedForFacility( i ) < scontext.fgraph.facilitiesPrices[i] - eps )
		{
			continue;
		}
		bool connFound = false;
		for ( int j = 0; j < facilitiesCount; j++ )
		{
			if ( scontext.specialSubgraph.existsEdge( i, j ) && scontext.indepFacilities[j] )
			{
				connFound = true;
				break;
			}
		}
		if ( !connFound )
		{
			result = false;
			buff << "\nvertex " << i << " is not in independent set and has no neighbors from the set";
		}
	}

	return result;
}

void TestComputeMaximalIndependentSet()
{
//	IGNORE( "TestComputeMaximalIndependentSet" );
	int n = 20;

	PrimDualSolvContext scontext( n, n );
	std::fill( scontext.fgraph.facilitiesPrices.begin(), scontext.fgraph.facilitiesPrices.end(), 0.0 );
	Graph::AdjMatrixType adjMatr = scontext.specialSubgraph.getAdjacentMatrixArr();

	// make a line
	for ( int i = 0; i < n; i++ )
	{
		int j = i ^ 1;
		adjMatr[i][j] = 1;
	}

	PrimalDualSolution::computeMaximalIndependentSet( scontext );

	BOOST_CHECK_EQUAL( 10, std::accumulate( scontext.indepFacilities.begin(), scontext.indepFacilities.end(), 0 ) );
	std::stringstream buff;
	bool independent = checkIndependent( scontext, buff );
	BOOST_CHECK( independent );
	if ( !independent )
	{
		std::cout << "\n" << buff.str();
	}

	// 20 .. 29 -- empty
	int n1 = 30;
	scontext.specialSubgraph.init( n1 );
	scontext.fgraph.facilitiesPrices.resize( n1, 0 );
	scontext.indepFacilities.resize( n1, false );
	adjMatr = scontext.specialSubgraph.getAdjacentMatrixArr();
	// make a line
	for ( int i = 0; i < n; i++ )
	{
		int j = i ^ 1;
		adjMatr[i][j] = 1;
	}
	// knows a lot
	int knowsAllNd2Vert = n;
	for ( int i = 0; i < n / 2; i++ )
	{
		adjMatr[knowsAllNd2Vert][i] = 1;
		adjMatr[i][knowsAllNd2Vert] = 1;
	}
	// extra edges in 10 .. 19
	for ( int i = n / 2; i < n; i++ )
	{
		if ( i % 3 != 1 ) continue;
		int j = i + 3;
		if ( j >= n ) continue;
		adjMatr[i][j] = 1;
		adjMatr[j][i] = 1;
	}
	std::fill( scontext.indepFacilities.begin(), scontext.indepFacilities.end(), false );

	PrimalDualSolution::computeMaximalIndependentSet( scontext );
	buff.clear();
	independent = checkIndependent( scontext, buff );
	BOOST_CHECK( independent );
	if ( !independent )
	{
		std::cout << "\n" << buff.str();
	}
	for ( int i = n + 1; i < n1; i++ )
	{
		BOOST_CHECK( scontext.indepFacilities[i] );
	}
}

void TestComputeMaximalIndependentSet_RandomCases()
{
//	IGNORE( "TestComputeMaximalIndependentSet_RandomCases" );
	int tcCount = BoostTestEnv::getMaximalIndependentSetTestCount();
	boost::variate_generator< boost::mt19937, boost::uniform_int<> > sampler(
		boost::mt19937( ( unsigned ) 12345 ), boost::uniform_int<>( 0, 1 << 30 ) );

	if ( tcCount > 5 ) std::cout << "\nTestComputeMaximalIndependentSet_RandomCases cases: " << tcCount;

	for ( int tc = 0; tc < tcCount; tc++ )
	{
		int n = 100;
		PrimDualSolvContext scontext( n, n );
		scontext.specialSubgraph.init( n );
		Graph::AdjMatrixType adjMatr = scontext.specialSubgraph.getAdjacentMatrixArr();

		for ( int i = 0; i < n; i++ )
		{
			for ( int j = 0; j < n; j++ )
			{
				int r = sampler() % 1000;
				if ( r < 200 && i != j )
				{
					adjMatr[i][j] = 1;
					adjMatr[j][i] = 1;
				}
			}
		}

		for ( int i = 0; i < n; i++ )
		{
			scontext.fgraph.facilitiesPrices[i] = 0;
			int r = sampler() % 1000;
			if ( r < 300 )
			{
				scontext.fgraph.facilitiesPrices[i] = 123;
				scontext.cityPays4fclt[i][i] = sampler() % 123;
			}
			if ( r < 150 )
			{
				scontext.cityPays4fclt[i][i] = 123;
				scontext.fgraph.setFSEdge( i, i, 0 );
			}
		}

		PrimalDualSolution::computeMaximalIndependentSet( scontext );
		std::stringstream buff;
		bool independent = checkIndependent( scontext, buff );
		BOOST_CHECK( independent );
		if ( !independent )
		{
			std::cout << "\n" << buff.str();
		}
	}
}

//TODO: order of increase cities payments is important !
void TestFacilityAssociation()
{
//	IGNORE( "TestFacilityAssociation" );
	try {
	FacilitiesGraph::FacilityPriceType EPS = 1e-5;
	int citiesCount = 6;
	int facilitiesCount = 4;
	PrimDualSolvContext scontext( citiesCount, facilitiesCount );

	scontext.fgraph.setFSEdge( 0, 0, 100 );
	scontext.fgraph.setFSEdge( 0, 1, 150 );
	scontext.fgraph.setFSEdge( 1, 0, 123 );
	scontext.fgraph.setFSEdge( 2, 0, 150 );
	scontext.fgraph.setFSEdge( 2, 1, 100 );
	scontext.fgraph.setFSEdge( 2, 2, 100 );
	scontext.fgraph.setFSEdge( 3, 1, 100 );
	scontext.fgraph.setFSEdge( 3, 2, 50 );
	scontext.fgraph.setFSEdge( 4, 2, 123 );
	scontext.fgraph.setFSEdge( 5, 3, 1234 );
	MakeSymmOperation::makeSymm( scontext.fgraph );

	scontext.fgraph.facilitiesPrices[0] = 100;
	scontext.fgraph.facilitiesPrices[1] = 150;
	scontext.fgraph.facilitiesPrices[2] = 150;
	scontext.fgraph.facilitiesPrices[3] = 4321;

	PrimalDualSolution::createWitnessForAllCities( scontext );
	BOOST_CHECK_EQUAL( 0, PrimalDualSolution::findWitnessForCity( scontext, 1 ) );
	BOOST_CHECK_EQUAL( 0, PrimalDualSolution::findWitnessForCity( scontext, 0 ) );
	BOOST_CHECK_EQUAL( 0, PrimalDualSolution::findWitnessForCity( scontext, 2 ) );
	BOOST_CHECK_EQUAL( 1, PrimalDualSolution::findWitnessForCity( scontext, 3 ) );
	BOOST_CHECK_EQUAL( 2, PrimalDualSolution::findWitnessForCity( scontext, 4 ) );
	BOOST_CHECK_EQUAL( 3, PrimalDualSolution::findWitnessForCity( scontext, 5 ) );

	PrimalDualSolution::computeSpecialEdgesSquareSubgraph( scontext );
	BOOST_CHECK_EQUAL( 4, scontext.specialSubgraph.getVertexCount() );
	BOOST_CHECK( scontext.specialSubgraph.existsEdge( 0, 1 ) );
	BOOST_CHECK( scontext.specialSubgraph.existsEdge( 2, 1 ) );
	BOOST_CHECK( !scontext.specialSubgraph.existsEdge( 2, 0 ) );

	BOOST_CHECK( !scontext.specialSubgraph.existsEdge( 0, 3 ) );
	BOOST_CHECK( !scontext.specialSubgraph.existsEdge( 1, 3 ) );
	BOOST_CHECK( !scontext.specialSubgraph.existsEdge( 2, 3 ) );

	PrimalDualSolution::computeMaximalIndependentSet( scontext );

	BOOST_CHECK_EQUAL( 3, std::accumulate( scontext.indepFacilities.begin(), scontext.indepFacilities.end(), 0 ) );
	BOOST_CHECK( scontext.indepFacilities[0] );
	BOOST_CHECK( scontext.indepFacilities[2] );
	BOOST_CHECK( scontext.indepFacilities[3] );

	BOOST_CHECK( PrimalDualSolution::isSpecialEdge( scontext, 0, 0 ) );

	PrimalDualSolution::computeFacilityAssociations( scontext );

	std::vector< FacilitiesGraph::FacilityPriceType > facilityEffPayment( facilitiesCount );
	for ( int i = 0; i < citiesCount; i++ )
	{
		int facility = scontext.associatedFacility[i];
		CRITICAL_CHECK( facility > -1 );
		BOOST_CHECK( scontext.indepFacilities[facility] );
		FacilitiesGraph::FacilityPriceType cityPays4Facility = scontext.cityPays4fclt[i][facility];
		facilityEffPayment[facility] += cityPays4Facility;
	}

	for ( int facility = 0; facility < facilitiesCount; facility++ )
	{
		if ( scontext.indepFacilities[facility] )
		{
			BOOST_CHECK_GE( facilityEffPayment[facility], scontext.fgraph.facilitiesPrices[facility] - EPS );
			if ( !( facilityEffPayment[facility] >= scontext.fgraph.facilitiesPrices[facility] - EPS ) )
			{
				std::cout << "\nfacility " << facility << " has cost: "
						<< scontext.fgraph.facilitiesPrices[facility] << " but payment was: "
						<< facilityEffPayment[facility];
			}
		}
	}

	} catch ( SolutionProcessException & exc )
	{
		std::cout << "\nTestFacilityAssociation, exception: " << exc.message;
	}
}

void TestFacilityAssociation_Random()
{
//	IGNORE( "TestFacilityAssociation_Random" );
	int tcCount = BoostTestEnv::getFacilityAssociationRandomTestCount();

	if ( tcCount > 5 ) std::cout << "\nTestFacilityAssociation_Random cases: " << tcCount;

	boost::variate_generator< boost::mt19937, boost::uniform_int<> > sampler(
		boost::mt19937( ( unsigned ) 12345 ),
		boost::uniform_int<>( 0, 1 << 30 )
	);

	CompoundGraphOperation * compOp = new CompoundGraphOperation();
	IGraphOperationPtr compOpPtr( compOp );

	GenBipartLinkNear * genNeigh = new GenBipartLinkNear();
	compOp->addSubOperation( IGraphOperationPtr( genNeigh ) );

	AddChpExpWeightsOperation * chpexp = new AddChpExpWeightsOperation();
	compOp->addSubOperation( IGraphOperationPtr( chpexp ) );

	compOp->addSubOperation( IGraphOperationPtr( new MakeSymmOperation() ) );

	GeneratePricesOperation * gpop = new GeneratePricesOperation();
	compOp->addSubOperation( IGraphOperationPtr( gpop ) );

	DefaultGraphGenerator ggen( compOpPtr );

	for ( int tc = 0; tc < tcCount; tc++ )
	{
		int citiesCount = 100 + sampler() % 100;
		int facilitiesCount = 30 + sampler() % 50;
		int neighbors = facilitiesCount / 10;
		neighbors += sampler() % neighbors;

		PrimDualSolvContext scontext( citiesCount, facilitiesCount );
		IOperationSettings oset( citiesCount, facilitiesCount );
		oset.neighborsBeg = neighbors;
		oset.neighborsRange = std::max( neighbors / 5, 1 );
		oset.cheapBeg = 100 + sampler() % 100;
		oset.cheapRange = 1 + sampler() % 20;
		oset.expensiveBeg = 300 + sampler() % 200;
		oset.expensiveRange = 1 + sampler() % 50;
		oset.facilityPriceBeg = 100 + sampler() % 300;
		oset.facilityPriceRange = 50 + sampler() % 100;
		ggen.generateGraph( scontext.fgraph, oset );
		BOOST_CHECK_EQUAL( citiesCount, scontext.fgraph.getCitiesCount() );
		BOOST_CHECK_EQUAL( facilitiesCount, scontext.fgraph.getFacilitiesCount() );

		bool created = PrimalDualSolution::createWitnessForAllCities( scontext );
		CRITICAL_CHECK( created );

		for ( int facility = 0; facility < facilitiesCount; facility++ )
		{
			FacilitiesGraph::FacilityPriceType price = scontext.fgraph.facilitiesPrices[facility];
			FacilitiesGraph::FacilityPriceType payed = scontext.getPayedForFacility( facility );
			BOOST_CHECK_LE( payed, price + 1e-5 );
		}

		for ( int city = 0; city < citiesCount; city++ )
		{
			int witness = PrimalDualSolution::findWitnessForCity( scontext, city );
			BOOST_CHECK_GT( witness, -1 );

			FacilitiesGraph::FacilityPriceType price = scontext.fgraph.facilitiesPrices[witness];
			FacilitiesGraph::FacilityPriceType payed = scontext.getPayedForFacility( witness );
			BOOST_CHECK_CLOSE( price, payed, 1e-3 );

			BOOST_CHECK_CLOSE(
				scontext.getCityPays4road( city, witness ),
				scontext.fgraph.getFSEdge( city, witness ),
				1e-3
			);


			for ( int facility = 0; facility < facilitiesCount; facility++ )
			{
				if ( scontext.fgraph.existsBiPartEdge( city, facility ) )
				{
					FacilitiesGraph::FacilityPriceType roadCost = scontext.fgraph.getFSEdge( city, facility );
					FacilitiesGraph::FacilityPriceType cityPays = scontext.getCityPays4road( city, facility );
					BOOST_CHECK_LE( cityPays, roadCost + 1e-5 );
				}
			}
		}
	}
}

void TestSolvContext__initFromGraph()
{
	PrimDualSolvContext scontext;
	BOOST_CHECK_EQUAL( 0, scontext.cityPays4fclt.size() );

	FacilitiesGraph gr( 5, 3 );
	scontext.fgraph = gr;
	BOOST_CHECK_EQUAL( 0, scontext.cityPays4fclt.size() );

	scontext.initFromGraph();
	BOOST_CHECK_EQUAL( 5, scontext.cityPays4fclt.size() );
	BOOST_CHECK_EQUAL( 3, scontext.cityPays4fclt[0].size() );
	BOOST_CHECK_EQUAL( 3, scontext.indepFacilities.size() );
	BOOST_CHECK_EQUAL( 5, scontext.associatedFacility.size() );
	BOOST_CHECK_EQUAL( 0.0, scontext.getCityPays4total( 5 - 1 ) );
}

void TestFactor3solver()
{
//	IGNORE( "TestFactor3solver" );
	FacilitiesGraph fgraph;

	IGraphGeneratorPtr graphGenerator(
		new DefaultGraphGenerator(
			CompoundGraphOperation::createOperation(
				4,
				new GenBipartLinkNear(),
				new AddChpExpWeightsOperation(),
				new MakeSymmOperation(),
				new GeneratePricesOperation()
			)
		)
	);
	DefaultGraphGenerator & stpg =
		dynamic_cast< DefaultGraphGenerator & >( * graphGenerator );
	CRITICAL_CHECK( & stpg != 0 );
	CompoundGraphOperation & stpop =
		dynamic_cast< CompoundGraphOperation & >( * stpg.getOperation() );
	CRITICAL_CHECK( & stpop != 0 );

	IFacilityLocationSolverPtr solver( new PrimalDualFacilitiesLocationSolver() );

	std::vector< char > isFacilityOpenVec;
	std::vector< int > cityWitnessVec;

	IOperationSettings oset( 50, 10 );
	oset.neighborsBeg = 5;
	oset.neighborsRange = 10;
	oset.facilityPriceBeg = 250;
	oset.facilityPriceRange = 300;
	oset.cheapBeg = 100;
	oset.cheapRange = 50;
	oset.expensiveBeg = 200;
	oset.expensiveRange = 100;
	graphGenerator->generateGraph( fgraph, oset );
	BOOST_CHECKPOINT( "pre-solve" );
	bool solved = solver->solve(
		fgraph,
		isFacilityOpenVec,
		cityWitnessVec
	);
	CRITICAL_CHECK( solved );
	CRITICAL_CHECK( isFacilityOpenVec.size() == fgraph.getFacilitiesCount() );
	CRITICAL_CHECK( cityWitnessVec.size() == fgraph.getCitiesCount() );

	BOOST_CHECKPOINT( "pre-check" );
	for ( int i = 0; i < fgraph.getCitiesCount(); i++ )
	{
		int witness = cityWitnessVec[i];
		BOOST_CHECK( witness > -1 );
		if ( witness == -1 ) break;
		BOOST_CHECK( fgraph.existsBiPartEdge( i, witness ) );
		BOOST_CHECK( witness < fgraph.getFacilitiesCount() );
		BOOST_CHECK( isFacilityOpenVec[witness] );
		if ( !isFacilityOpenVec[witness] ) break;
	}
	FacilitiesGraph::FacilityPriceType totalPDCost = IFacilityLocationSolver::computeCost(
		fgraph,
		isFacilityOpenVec,
		cityWitnessVec
	);

	std::vector< char > isFacilityOpenVecB;
	std::vector< int > cityWitnessVecB;
	BinaryRunFacilityLocationSolver().solve(
		fgraph,
		isFacilityOpenVecB,
		cityWitnessVecB
	);
	FacilitiesGraph::FacilityPriceType totalPDCostB = IFacilityLocationSolver::computeCost(
		fgraph,
		isFacilityOpenVecB,
		cityWitnessVecB
	);

	BOOST_CHECK_LT( totalPDCost, totalPDCostB * 3 );

	// ### add triangle
	stpop.addSubOperation( IGraphOperationPtr( new TriangleInequalityCloseOperation() ) );

	oset.firstBipart = 100;
	oset.secondBipart = 15;
	oset.neighborsBeg = 5;
	oset.neighborsRange = 10;
	graphGenerator->generateGraph( fgraph, oset );
	solved = solver->solve(
		fgraph,
		isFacilityOpenVec,
		cityWitnessVec
	);
	totalPDCost = IFacilityLocationSolver::computeCost(
		fgraph,
		isFacilityOpenVec,
		cityWitnessVec
	);
	BinaryRunFacilityLocationSolver().solve(
		fgraph,
		isFacilityOpenVecB,
		cityWitnessVecB
	);
	totalPDCostB = IFacilityLocationSolver::computeCost(
		fgraph,
		isFacilityOpenVecB,
		cityWitnessVecB
	);
	BOOST_CHECK_LT( totalPDCost, totalPDCostB * 2 );
}

void TestFactor3solver_Random()
{
//	IGNORE( "TestFactor3solver_Random" );
	int smallTestCaseCount = BoostTestEnv::getFactor3RandomSmallTestCount();
	int testCaseCount = BoostTestEnv::getFactor3RandomTestCount();
	bool printCostReport = BoostTestEnv::getFactor3RandomPrintCost();

	std::stringstream devnull;
	std::ostream & costsReport = printCostReport ? std::cout : devnull;

	if ( testCaseCount > 100 ) std::cout << "TestFactor3solver_Random test cases: " << testCaseCount;

	boost::variate_generator< boost::mt19937, boost::uniform_int<> > sampler(
		boost::mt19937( ( unsigned ) 12345 ),
		boost::uniform_int<>( 0, 1 << 30 )
	);

	FacilitiesGraph fgraph;

	IGraphGeneratorPtr graphGenerator(
		new DefaultGraphGenerator(
			CompoundGraphOperation::createOperation(
				5,
				new GenBipartLinkNear(),
				new AddChpExpWeightsOperation(),
				new MakeSymmOperation(),
				new GeneratePricesOperation(),
				new TriangleInequalityCloseOperation()
			)
		)
	);
	DefaultGraphGenerator & stpg =
		dynamic_cast< DefaultGraphGenerator & >( * graphGenerator );
	CRITICAL_CHECK( & stpg != 0 );
	CompoundGraphOperation & stpop =
		dynamic_cast< CompoundGraphOperation & >( * stpg.getOperation() );
	CRITICAL_CHECK( & stpop != 0 );

	IFacilityLocationSolverPtr solver( new PrimalDualFacilitiesLocationSolver() );
	std::vector< char > isFacilityOpenVec;
	std::vector< int > cityWitnessVec;
	std::vector< char > isFacilityOpenVecB;
	std::vector< int > cityWitnessVecB;

	for ( int tc = 0; tc < testCaseCount; tc++ )
	{
		if ( tc % 3 == 0 ) costsReport << "\n"; else costsReport << "\t";

		bool isSmallTC = tc < smallTestCaseCount;

		int citiesCount = -1;
		int facilitiesCount = -1;
		if ( isSmallTC )
		{
			citiesCount = 50 + sampler() % 50;
			facilitiesCount = 10 + sampler() % 3;
		} else
		{
			citiesCount = 300 + sampler() % 150;
			facilitiesCount = 100 + sampler() % 50;
		}

		IOperationSettings oset( citiesCount, facilitiesCount );
		oset.neighborsBeg = 5;
		oset.neighborsRange = 10;
		oset.cheapBeg = 100;
		oset.cheapRange = 50;
		oset.expensiveBeg = 200;
		oset.expensiveRange = 100;
		oset.facilityPriceBeg = 200;
		oset.facilityPriceRange = 100;
		graphGenerator->generateGraph( fgraph, oset );

		bool solved = solver->solve(
			fgraph,
			isFacilityOpenVec,
			cityWitnessVec
		);
		CRITICAL_CHECK( solved );
		CRITICAL_CHECK( isFacilityOpenVec.size() == fgraph.getFacilitiesCount() );
		CRITICAL_CHECK( cityWitnessVec.size() == fgraph.getCitiesCount() );

		for ( int i = 0; i < fgraph.getCitiesCount(); i++ )
		{
			int witness = cityWitnessVec[i];
			CRITICAL_CHECK( witness > -1 );
			BOOST_CHECK( fgraph.existsBiPartEdge( i, witness ) );
			CRITICAL_CHECK( witness < fgraph.getFacilitiesCount() );
			BOOST_CHECK( isFacilityOpenVec[witness] );
			if ( !isFacilityOpenVec[witness] ) break;
		}
		FacilitiesGraph::FacilityPriceType totalPDCost = IFacilityLocationSolver::computeCost(
			fgraph,
			isFacilityOpenVec,
			cityWitnessVec
		);
		costsReport << "tc: " << tc << " [f3r = " << totalPDCost << " best = ";
		if ( isSmallTC )
		{
			BinaryRunFacilityLocationSolver().solve(
				fgraph,
				isFacilityOpenVecB,
				cityWitnessVecB
			);
			FacilitiesGraph::FacilityPriceType totalPDCostB = IFacilityLocationSolver::computeCost(
				fgraph,
				isFacilityOpenVecB,
				cityWitnessVecB
			);
			BOOST_CHECK_LT( totalPDCost, totalPDCostB * 3 );
			costsReport << totalPDCostB;
		} else
		{
			costsReport << "?";
		}
		costsReport << "]";
	}
}
